/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#ifndef __SQUARED_LOSS_HH__
#define __SQUARED_LOSS_HH__

#include <boost/optional.hpp>
#include <cmath>
#include <stdint.h>
#include "loss.hh"

namespace flassol
{
  class SquaredLoss : public Loss
    {
      private:
        static LossFactory factory;

        double   sumy;
        double   sumysq;
        uint64_t n;

      public:
        explicit
        SquaredLoss () : sumy (0), sumysq (0), n (0) { }

        virtual float 
        loss (float p,
              float y)
          {
            sumy += y;
            sumysq += y * y;
            ++n;

            float z = (p - y);
            return 0.5 * z * z;
          }

        virtual float 
        dldp (float p,
              float y)
          {
            return (p - y);
          }

        virtual float 
        invariant_update (float p,
                          float wdotx,
                          float adotb,
                          float y,
                          float eta,
                          float anormplusbnorm,
                          float xnorm,
                          float onedotone)
          {
            // TODO: handle dyadic case
            (void) wdotx;
            (void) adotb;
            (void) anormplusbnorm;
            (void) onedotone;

            // The actual learning rate is (eta / norm)
            // So exponent is (norm * eta / norm) = eta
            
            if (fabs (eta) < 1e-5)
              {
                return (p - y) * eta / xnorm;
              }
            else
              {
                return (p - y) * (1.0f - expf (-eta)) / xnorm;
              }
          }

        virtual std::ostream& 
        report (std::ostream& o,
                LogLevel      level)
          {
            if (n > 0 && level <= INFO)
              {
                double var = (sumysq / n) - (sumy / n) * (sumy / n);

                o << "best constant = " << (sumy / n) << std::endl;
                o << "best constant loss = " << 0.5 * var << std::endl;
              }

            return o;
          }
    };

  class SquaredClippedLoss : public Loss
    {
      private:
        static LossFactory factory;

        SquaredLoss squared_loss;
        boost::optional<float> minlabel;
        boost::optional<float> maxlabel;

        void
        update_bounds (float y)
          {
            if (! maxlabel || y > *maxlabel)
              {
                maxlabel = y;
              }

            if (! minlabel || y < *minlabel)
              {
                minlabel = y;
              }
          }

        float
        clip (float p)
          {
            return (p > *maxlabel) ? *maxlabel :
                   (p < *minlabel) ? *minlabel : p;
          }

      public:
        explicit
        SquaredClippedLoss () { }

        explicit
        SquaredClippedLoss (float min, 
                            float max) 
          : minlabel (min),
            maxlabel (max)
          {
          }

        virtual float loss (float p,
                            float y)
          {
            update_bounds (y);
            return squared_loss.loss (clip (p), y);
          }

        virtual float dldp (float p,
                            float y)
          {
            update_bounds (y);
            return squared_loss.dldp (clip (p), y);
          }

        virtual float 
        invariant_update (float p,
                          float wdotx,
                          float adotb,
                          float y,
                          float eta,
                          float anormplusbnorm,
                          float xnorm,
                          float onedotone)
          {
            update_bounds (y);
            return squared_loss.invariant_update (clip (p), 
                                                  wdotx,
                                                  adotb,
                                                  y, 
                                                  eta, 
                                                  anormplusbnorm,
                                                  xnorm,
                                                  onedotone);
          }

        virtual std::ostream& 
        report (std::ostream& o,
                LogLevel      level)
          {
            if (minlabel && level <= INFO)
              {
                o << "minimum label encountered = " << *minlabel << std::endl;
              }

            if (maxlabel && level <= INFO)
              {
                o << "maximum label encountered = " << *maxlabel << std::endl;
              }

            return squared_loss.report (o, level);
          }
    };
}

#endif // __SQUARED_LOSS_HH__
